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We investigate the quantum phase transition in an S = 1/2 dimerized Heisenberg antiferromagnet 
in three spatial dimensions. By performing large-scale quantum Monte Carlo simulations and de¬ 
tailed finite-size scaling analyses, we obtain high-precision results for the quantum critical properties 
at the transition from the magnetically disordered dimer-singlet phase to the antiferromagnetically 
ordered Neel phase. This transition breaks O(IV) symmetry with A = 3inD = 3-|-l dimensions. 

This is the upper critical dimension, where multiplicative logarithmic corrections to the leading 
mean-field critical properties are expected; we extract these corrections, establishing their precise 
forms for both the zero-temperature staggered magnetization, m s , and the Neel temperature, Tjv- 
We present a scaling Ansatz for Tjv, including logarithmic corrections, which agrees with our data 
and indicates exact linearity with m 3 , implying a complete decoupling of quantum and thermal 
fluctuation effects even arbitrarily close to the quantum critical point. We also demonstrate the 
predicted A-independent leading and subleading logarithmic corrections in the size-dependence of 
the staggered magnetic susceptibility. These logarithmic scaling forms have not previously been 
identified or verified by unbiased numerical methods and we discuss their relevance to experimental 
studies of dimerized quantum antiferromagnets such as TlCuCB. 

PACS numbers: 75.10.Jm, 75.40.Cx, 75.40.Mg 


I. INTRODUCTION 

Antiferromagnetic insulators exhibit a multitude of 
fundamental phenomena in the neighborhood of the 
phase transitions separating their magnetically ordered 
ground states from different types of quantum paramag¬ 
netic phase. These quantum phase transitions (QPTs) 
occur at temperature T = 0 as a consequence of non- 
thermal parameters (examples include magnetic fields, 
applied pressure, and dopant concentration) that act to 
change the effect of quantum mechanical fluctuations 
mm- At finite temperatures, a further dimension is 
opened in the presence of both quantum and classical 
(thermal) fluctuations, and the rich physics arising from 
their interplay includes all the properties of the quantum 
critical (QC) regime [3]. 

Experimentally, the material in which the most de¬ 
tailed study of intertwined classical and quantum criti¬ 
cal behavior has been performed is TICUCI 3 . This com¬ 
pound is composed of antiferromagnetically coupled pairs 
of Cu 2+ ions (S = 1/2), which tend to form dimer sin¬ 
glets and have antiferromagnetic interdimer couplings in 
all three spatial dimensions (d = 3) [4]. At ambient pres¬ 
sure and zero field, TICUCI 3 is a nonmagnetic insulator 
with a gap of 0.63 meV to triplet spin excitations. As a 
consequence of this small gap, an applied magnetic field 
of 5.4 T is sufficient to drive the system to an ordered 
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antiferromagnetic state, through a QPT in the Bose- 
Einstein universality class j5]. A relatively small applied 
hydrostatic pressure, p c = 1.07 kbar, is also sufficient to 
create an antiferromagnetically ordered state [B], through 
a QPT in the three-dimensional (3D) 0(3) universality 
class due to spontaneous breaking of the SU(2) spin sym¬ 
metry (which is further reduced in TICUCI 3 by a weak 
unixial anisotropy, making the universality class 3D XY). 

The elementary excitations on the ordered side of the 
zero-field quantum critical point (QCP) are gapless spin 
waves, the Goldstone modes associated with spontaneous 
breaking of spin-rotational symmetry. On the disordered 
side, quantum fluctuations, towards spin-singlet forma¬ 
tion on the dimers, suppress the long-range antiferromag¬ 
netic order, restoring the symmetry and ensuring that 
all excitations are gapped. This evolution of the excita¬ 
tion spectrum in TIC 11 CI 3 has been measured in Ref. [7] . 
At finite temperatures on the ordered side, a classical 
phase transition occurs at the Neel temperature, Tjv, 
where the long-ranged magnetic order is “melted” not 
by quantum but by thermal fluctuations. At finite tem¬ 
peratures around the QCP, the combination of quantum 
and thermal fluctuations creates the QC regime, where 
the only characteristic energy scale of the system is the 
temperature itself and many universal properties emerge 
[3]. The phase diagram of TICUCI 3 under pressure and 
the restoration of classical critical scaling around X)v were 
the subject of a recent investigation [5 . 

QPTs in dimerized quantum spin models have been 
studied numerically by a number of authors, primarily 
by quantum Monte Carlo (QMC) simulations. Early in- 
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FIG. 1. (Color online) (a) Dimerized lattice of S = 1/2 spins 
in the 3D double cubic geometry. Sites of the red and blue cu¬ 
bic lattices are connected pairwise by dimer bonds; J' and J 
are antiferromagnetic Heisenberg interactions, respectively on 
and between the dimer units, and their ratio, g = J'/J, con¬ 
trols the QPT from a Neel ordered phase (left) to a quantum 
disordered dimer-singlet phase (right), with the QCP occur¬ 
ring at the critical ratio g c . (b) Schematic quantum critical 
phase diagram for the Heisenberg model on the double cubic 
lattice. On the ordered “renormalized classical” side, g < g c , 
the Neel order is progressively weakened by increasing quan¬ 
tum fluctuations as g approaches g c , causing both the order¬ 
ing (Neel) temperature, Tjv(g), and the order parameter (the 
staggered magnetization), m s (g), to go continuously to zero. 
On the “quantum disordered” side, intradimer correlations 
dominate and the characteristic energy scale, A, is the gap to 
triplet excitations. Above g c at T > 0 is the universal “quan¬ 
tum critical” (QC) region, whose behavior is governed by the 
3 + 1-dimensional 0(3) universality class. Around Tjv(<7), one 
expects a region of classical critical (CC) behavior where ther¬ 
mal fluctuations are dominant. 


vestigations of the bilayer square-lattice antiferromagnet 
[5] and other two-dimensional geometries [TU] HI] have 
been followed more recently by high-precision studies of 
a range of critical properties [T2MT6] . In three dimensions, 
the focus of investigations has been on the field-induced 
transition na. on the effects of dimensionality iisiiis], 
and on physical observables at the coupling-induced QPT 

PH- 

A minimal model of a dimerized quantum antiferro¬ 
magnet has only two coupling constants, J' on and J 


between the dimer units, and therefore only one con¬ 
trol parameter, g = J'/J. The geometry considered in 
the present study is the double cubic lattice shown in 
Fig. 0 a). In this system at large g , intradimer singlet 
correlations dominate the physics and the ground state 
is magnetically disordered, while at small g the inter- 
dimer correlations establish long-ranged magnetic order. 
The order parameter of the Neel phase is the staggered 
magnetization, m s (g), and along with the ordering tem¬ 
perature, Tjv(g), it can be driven continuously to zero by 
increasing g 1 as illustrated in Fig. 0 b). By standard ar¬ 
guments of dimensionality and symmetry, the dynamical 
exponent of this system is 2 = 1 0 [23] and the QCP 
belongs to the D = 3 + 1 0(3) universality class, which 
is at the upper critical dimension, D c = 4, of all O (N) 
models Gann. At D = D c , mean-field critical scaling 
behavior alone is not sufficient to capture the physics 
of fluctuations around the QCP and multiplicative log¬ 
arithmic corrections to the physical quantities (thermo¬ 
dynamic functions) are expected. 

The theoretical importance of multiplicative logarith¬ 
mic corrections to mean-field scaling behavior lies not 
only in the statistical physics of condensed matter sys¬ 
tems but also in high-energy physics [MIHS]- The general 
problem of a quantum field theory with an iV-component 
field is encapsulated by an “O(iV) </> 4 theory,” a La- 
grangian containing a dynamic (quadratic gradient) term 
and a potential term with quadratic ( (f > 2 ) and quartic ((/ 4 ) 
contributions. On changing the sign of the quadratic 
term, the system is driven through a QPT separating a 
phase with (<f>) = 0 from one with ( (f >) 7 ^ 0 (a “Mexican 
hat” potential). As noted above, the low-energy proper¬ 
ties of the 3D dimerized antiferromagnet of SU(2) quan¬ 
tum spins with Heisenberg interactions correspond to a 
field theory with N = 3 and D = 3 + 1 (including the 
time dimension); N = 1 and 2 correspond respectively to 
Ising and XY spin interactions. 

Beyond the upper critical dimension ( D > D c ), the 
scaling behavior of the 0(N) </> 4 theory is straightfor¬ 
ward, with the critical exponents being exactly those 
given by mean-field theory pSl - ET] . namely a = 0 , /3 = 
1/2, 7 = 1 , 6 = 3, and v = 1/2. As we discuss below, 
this may be taken as an expression of the independence of 
quantum and thermal fluctuations when the phase space 
is sufficiently large. For D < D c , the situation is complex 
and these exponents take anomalous values. However, 
exactly at the upper critical dimension, D = D c = 4, the 
leading scaling behavior coincides with that of the mean- 
field theory, but modified by multiplicative logarithmic 
corrections [211 j2SJ EH] ■ While the leading exponents are 
IV-independent, a measure of iV-dependence resides in 
the logarithmic corrections and these must be taken into 
account to establish the universality class of the transi¬ 
tion [SO]. Because the established results for the form 
of these corrections are based on perturbative techniques 
applied to low-energy theories, it is desirable to verify 
them using unbiased numerical methods applied directly 
to the lattice Hamiltonians, and this is what we achieve 
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here. 

Despite the insight into general QPT phenomena ob¬ 
tained from simulations using this type of minimal model 
for 3D dimerized systems the question of 

logarthmic corrections has to date been addressed only 
briefly and inconclusively mm- Experimentally, the 
feasibility of observing logarithmic corrections in systems 
such as TICUCI 3 remains a challenging open issue [5] . In 
this paper we provide a systematic numerical study. We 
employ large-scale QMC simulations to investigate the 
critical behavior of the order parameter and Neel tem¬ 
perature on the double cubic lattice [Fig. [l|a)] for small 
values of \g — g c \ unattainable in all previous studies. 
State-of-the-art QMC techniques and finite-size-scaling 
analysis using very large systems (sizes exceeding 10 5 
spins) allow us to detect and characterize the multiplica¬ 
tive logarithmic corrections in the universal scaling rela¬ 
tions for the QC regime at the upper critical dimension, 
here for the D = 3 + 1 0(3) QCP. In fact our results con¬ 
stitute hitherto unavailable exact numerical verification 
of the logarithmic forms predicted both by perturbative 
renormalization-group calculations based on the 

iV-component <(> 4 theory at D c and by additional con¬ 
siderations exploiting the zeros of the partition function 
[25l l28l l29l 135] . To the best of our knowledge, no system¬ 
atic numerical calculations have been performed beyond 
N = 1 [ 251155 ]. 

As will become clear, our numerical results demon¬ 
strate to high precision the validity of the detailed the¬ 
oretical predictions for the expected universality class. 
Both size-dependent scaling and the order parameter in 
the thermodynamic limit show evident deviations from 
pure mean-field behavior, which are accounted for by log¬ 
arithmic corrections whose exponents are in very close 
ageement with the predicted values where available. In 
the case of the Neel temperature, we are not aware of any 
previous scaling predictions including logarithms. Here 
we test an Ansatz based on the known scaling behavior 
for the relevant energy scales [3133 and the logarithmic 
corrections in corresponding classical systems [25] ■ In 
addition to probing the asymptotic behavior, our results 
also provide direct insight into the range of validity of log¬ 
arithmically modified critical scaling forms as one moves 
away from the QCP, which will be essential in evaluating 
the experimental relevance of logarithmic corrections. 

The paper is organized as follows. In Sec. [H] we intro¬ 
duce the model and the numerical method, describing the 
measurement of physical observables in our QMC sim¬ 
ulations. In Sec. |III| we begin the presentation of our 
numerical results with the precise determination of g c , 
the position of the QCP, using finite-size-scaling tech¬ 
niques. Section |TV| discusses the observation of clear log¬ 
arithmic corrections in the staggered magnetic suscepti¬ 
bility, x(Qafj L), at the QCP as a function of the sys¬ 
tem size, L. We present our results for the sublattice 
magnetization, m s , at T = 0 in Sec.[Vj discussing in de¬ 
tail its extrapolation to the thermodynamic limit, where 
we investigate the presence of logarithmic corrections to 


the leading mean-field behavior. In Sec. |VI| we present a 
scaling Ansatz for the Neel temperature, apply finite-size 
scaling to extract it as a function of g 1 and again inves¬ 
tigate corrections to mean-field behavior. We compute 
the characteristic velocity c of spin excitations, demon¬ 
strate the precise linearity of T/v and m s , and discuss 
the physical interpretation of this behavior. We summa¬ 
rize our results in Sec. IVIll and comment further on their 
theoretical and experimental consequences. 


II. MODEL AND METHODS 

As a representative 3D dimerized lattice with an un¬ 
frustrated geometry, we choose to study the double cu¬ 
bic model shown in Fig. 0 a )- This system consists of 
two interpenetrating cubic lattices with the same anti¬ 
ferromagnetic interaction strength, J, connected pairwise 
by another antiferromagnetic interaction, J'. The QPT 
occurs when the coupling ratio g = J' / J is increased, 
changing the ground state from a Neel-ordered phase of 
finite staggered magnetization to a dimer-singlet (“quan¬ 
tum disordered”) phase, as illustrated in Fig. 0b). An 
advantage of this geometry over cases where the dimer¬ 
ization is imposed within a single lattice, such as the 
simple cubic lattice m, is that all symmetries of the cu¬ 
bic lattices are retained, facilitating the consideration of 
quantities such as the spin stiffness or the velocity of spin 
excitations. 

The Hamiltonian is given by 

h = y ' JijSi ■ §j, ( 1 ) 

(ij) 

where S) is an S = 1 /2 spin operator residing on a double 
cubic lattice of N = 2 L 3 sites with periodic boundary 
conditions. The sum is taken only over nearest-neighbor 
sites, where every site has six neighbors on the same cubic 
lattice with coupling strength = J and one neighbor 
on the opposite cubic lattice with J, ? = J' [Fig. JlTa)] . 
We set J = 1 as the unit of energy and use g = Jy J as 
the control parameter. 

To study this system, we use the stochastic series 
expansion (SSE) QMC technique fl5j 138H40] to ob¬ 
tain unbiased results, i.e. numerically exact within well- 
characterized statistical errors, for physical quantities in 
systems of finite side-length L. Here we present results up 
to L = 48 at temperatures T = j3 _1 with /3 up to 2 L. We 
then perform detailed analyses by finite-size scaling [41 j 
to extract information in the thermodynamic limit both 
in the ordered state and at the QCP, as detailed in the 
separate sections to follow. Here we define the physical 
quantities of interest and discuss some technical aspects 
of their calculation within the SSE QMC method. 

Because spin-rotation symmetry is not broken in sim¬ 
ulations of finite-size systems, one may measure the 
squared order parameter and take its square root as a 
post-simulation step. The staggered magnetization is 
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given by 

m s = \J ^S{ Qaf), (2) 


which is dimensionless and satisfies the crucial property 
of being size-independent at the QCP in the limit of large 
system sizes. The spin stiffness, or helicity modulus, of 
the system is defined as 


where 

1 N 

S(<l) = ^e- i ^-^(S i -S j ) (3) 

hi 

is the magnetic structure factor, with denoting the 
real-space position of the spin Si on lattice site i, and 
Qaf = (tt, 7t, 7t, n) is the wave vector of antiferromag¬ 
netic order, with the fourth 7r denoting the phase factor 
between the two cubic lattices. We consider only the , 0 - 
component of the magnetization and average it over the 
time dimension of the QMC configurations, computing 
the expectation value of the squared order parameter in 
the form 


Ps = 


1 d 2 F^ 0 


N d 2 (j) 0 


a = x,y,z, 


( 10 ) 


where F is the free energy and (f> a is the angle of a twist 
imposed between all spins in planes perpendicular to the 
a axis. In an SSE simulation, the most efficient way 
to extract the spin stiffness is to take the derivative in 


Eq. (10) directly in the QMC expression for F((j> a ) at 


= 0, giving 


P 


a. 

s 


3(Wq) 

4/3 


( 11 ) 


where 


m 


2 

sz 


1 

p 



dTm 2 sz (r ), 


(4) 


Wa = \{N+ - N-) 


( 12 ) 


where 

™ sz (T) = if>- iQAF ' ri Sf(T), (5) 

V 2=1 

with 

Si( T ) = e rH S^e~ rH (6) 

the time-evolved spin operator at imaginary time r. In an 
SSE simulation, the integral in Eq. Q is transformed into 
a discrete sum with no approximations and the relation 
compensating for the rotational averaging of the single 
measured component of the order parameter, 

m s = (7) 

is applied post-simulation. The magnetic susceptibility 
is defined as 

X(q) = ^11^ dT(S?(T)S*(0))e-^-^. (8) 


is the winding number [42l [44] of the spin in spatial direc¬ 
tion a and 7V+ and N~ are the numbers of occurrences 
of the operators SP SJ and S~ S+ on bond (i, j) in the a- 
direction within imaginary time [0,/3]. As noted above, 
is the same in all three directions due to the sym¬ 
metry of the double cubic lattice, and the average may 
be taken over all of these. The spin stiffness follows the 
scaling law p s oc L 2 ~ d ~ z in d spatial dimensions pA] and, 
because the dynamic exponent is 2 = 1 here, the quantity 
p s L d ~ , or equivalently p s L D ~ 2 , is also size-independent 
at the QCP, up to a logarithmic correction at the upper 
critical dimension. 

Finally, the spin-wave velocity, c, can be obtained re¬ 
liably by monitoring the fluctuations of the spatial and 
temporal winding numbers [E1I221030B]. For a fixed 
system size, L, the inverse temperature /3 is adjusted to 
the value /3* where the system has equal winding-number 
fluctuations in the spatial and temporal directions, 

(M’q(/3*)) = ( w t(P*))i a = x,y, z. (13) 


In the SSE approach, the squared order parameter 
[Eq. Q] is readily evaluated at any r because the QMC 
configurations are constructed in the S z basis EH 02]. 
The dynamical spin-spin correlation function contained 
in Eq. ([8]) can also be obtained easily by applying an op¬ 
erator string connecting the S z states at different imag¬ 
inary times, with the integral computed analytically to 
give a direct, formally exact QMC estimator not requir¬ 
ing post-simulation integration mM- 

The Binder ratio [43] is the ratio of the fourth moment 
of a quantity to the square of its second moment. For 
our purposes, the relevant quantity is 


R, 


{miz) 

(™ 2 Z ) 2 ’ 


(9) 


The temporal winding number is the net magnetization 
(number of up minus number of down spins) of the sys¬ 
tem, w T = M z = S z , which is easily obtained in the 
S z basis [16j- When the condition is met, the spin- 
wave velocity is given by the ratio 


L 

Ml)' 


(14) 


The isotropy of the lattice is an advantage also in this 
case. For each value of g , an extrapolation L —> 00 is 
performed to obtain c in the thermodynamic limit. For 
further details of these procedures, we refer the reader 
to the recent extensive tests of this method conducted in 
Ref. [161. 
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III. DETERMINATION OF THE QCP 

The key to an accurate characterization of logarithmic 
corrections is a high-precision determination of the loca¬ 
tion g c of the QCP. For this we employ the Binder ratio, 
i? 2 i and the appropriately scaled spin stiffness, p s L D ~ 2 , 
which both have scaling dimension zero and therefore 
should approach constant values at g c when L —> oo, up 
to possible logarithmic corrections. We stress that the 
scaling forms for the approach of both quantities to the 
critical point are valid for a four-dimensional (4D) the¬ 
ory, with the temperature (imaginary time) providing the 
fourth dimension, and are applicable on a “critical con¬ 
tour” where the inverse temperature T -1 = kL is taken 
to infinity symmetrically with the spatial dimension of 
the system. This form is appropriate for a system with 
dynamic exponent z =1 (in general 1/T ~ L z ). All val¬ 
ues of k yield the same results in the limit L —> oo, and in 
principle the contour is optimal when k = 1 /c; the spin- 
wave velocity, c, is a number of order unity discussed in 
detail in Sec. m and here we use k = 1 . 

Away from g c , R 2 and p s L 2 approach different con¬ 
stant values with increasing system size. In the Neel 
state, R 2 approaches 9/5 due to diminishing fluctuations 
in the magnitude of the rotationally-invariant order pa¬ 
rameter, of which we measure only the ^-component in 
Eq.@. In the quantum disordered phase, R 2 approaches 
a higher value dictated by Gaussian fluctuations, which 
from the symmetries of the double cubic model is 3. The 
spin stiffness falls from non-zero values in the Neel phase 
to zero in the disordered phase. When calculated as func¬ 
tions of g , the curves R 2 (g,L) and p s (g,L)L 2 obtained 
for different system sizes should cross at the QCP, up to 
corrections that are well understood from the theory of 
finite-size scaling. We analyze these corrections to ob¬ 
tain an unbiased value of the critical coupling, g Cl in the 
thermodynamic limit m- 

Figure [2](a) shows R 2 as a function of g in the neigh¬ 
borhood of the critical coupling ratio for various system 
sizes. We have performed simulations for all even-length 
sizes L = 6, 8 , 10, ..., 40, but here we present only the 
L = 30, 32, ..., 40 data for clarity. Analogous curves for 
the scaled spin stiffness, p s L 2 , are shown in Fig. Hb), 
again only for system sizes L = 30, 32, ..., 40. In 
both cases, the system sizes are sufficiently large that 
the crossing points exhibit only a very weak dependence 
on L on the scale used in the figure, and both sets of data 
may be used independently to show that the QCP is lo¬ 
cated at g c ~ 4.837(1). A detailed analysis is required 
to obtain the most precise results attainable, free of any 
finite-size effects, and we first discuss the general scaling 
behavior of g c before presenting our numerical results. 

A. Scaling Forms for Critical-Point Estimators 

To describe the evolution of the crossing points with 
L, we perform a systematic extrapolation of the finite- 




1/L 


FIG. 2. (Color online) (a) Binder ratio, R 2 , and (b) scaled 
spin stiffness, p s L 2 , as functions of the coupling ratio g for 
system sizes L = 30, 32, ..., 40. Crossings of the curves for 
pairs of system sizes L and L + 2 define finite-size estimates 
g^(L) and g£(L) of the critical point, which are are fitted 
to the form of Eq. ( |16| ) in panel (c). Enforcing a common 
value of g c in both fits gives the L —> 00 critical point as 
g c = 4.83704(6) and the irrelevant exponent as lo = —0.31(5) 
for i ?2 and u> = 0.82(5) for p s L 2 . 


size data to the thermodynamic limit by extracting the 
crossing points between data sets for all pairs of system 
sizes, L and L + 2, based on polynomial interpolations. 
Figure [2](c) shows the crossing points <?,?(£) and g£(L) 
obtained in this manner. 

For any quantity probing a singularity in the thermo¬ 
dynamic limit, one may define a size-dependent critical 
point g' c (L). In general, this quantity is expected to shift 
by an amount proportional to L -1 /" with respect to the 
true infinite-size critical point, g c (p o) (hereafter denoted 
for simplicity by g c ), i.e. for large L , 

9c(L) = g c + aL ~ 1/v , (15) 

where v is the standard correlation-length exponent. 
However, with a definition based on crossing points of a 
dimensionless quantity computed for two different sizes, 
the leading corrections cancel and the convergence is 
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faster, 


g c (L)=g c + aL~^ 1 '+ u \ 


(16) 


where to > 0 is the dominant irrelevant exponent. In 
practice, with data fits to a rather limited range of avail¬ 
able system sizes, the corrections to Eq. (15) contained in 


Eq. (16) will have exponents and prefactors that deviate 


from their asymptotic values due to the neglected correc¬ 
tions of higher order in 1/L and therefore these should 
be considered as “effective” quantities. 

The above forms are applicable in the absence of loga¬ 
rithmic corrections, but such corrections are the primary 
focus of our study and are expected at D = D c . Kenna 
has derived the modified form of Eq. (151 for classical 
systems with logarithmic corrections [251 147j . 


g c (L) = g c + aL 1,v ln A L, 


(17) 


neglected higher-order irrelevant fields, and the expected 
weak logarithmic corrections. However, the good match 
obtained between the two extrapolated <? c -estimators, es¬ 
pecially when approaching the infinite-size value from dif¬ 
ferent directions, would not be expected in the presence 
of any corrections not taken sufficiently into account by 
the fitting functions. Thus we believe the error bar on 
g c quoted above to be completely representative of all 
statistical and systematic uncertainties, in the sense that 
any remaining systematic errors due to the fitting form 
should be smaller than the statistical errors. The tests we 
perform on the critical scaling behavior around the QCP 
in the subsequent sections also support this statement. 


IV. SIZE-DEPENDENT LOGARITHMIC 
CORRECTIONS AT THE QCP 


where the exponent of the logarithm for the 4D 0(3) 
universality class is A = —1/22. For the crossing points, 
Eq. (16) is modified to 


g c (L) = g c + aL-^ v +^ ln £ L, 


(18) 


as shown in Appendix A, where c = A if the subleading 
term L ~“ has no multiplicative logarithmic correction, 
but is altered by an unknown amount if it does. Un¬ 
der the circumstances, with a number of unknowns and 
with simulation data only for a restricted range of system 
sizes, we fit our data using not Eq. (18 1 but instead the 


purely algebraic form of Eq. (16 1 with v = 1/2 and w, 
the effective value of the subleading exponent over the 
fitting range, treated as a different fitting parameter for 
the separate quantities i ?2 and p s L 2 . 


B. Numerical Determination of g c 

We take both pairs of fitting parameters a and ui in 
Eq. ( |T6| ) to be independently free for the two datasets 
g^{L)~ and g£{L), but impose the constraint that the 
curves have the same g c . As shown in Figs. [ 2 ]) a) and 
DM, we obtain good fits to both functions, meaning with 
a reduced y 2 (per degree of freedom, hereafter denoted 
Xr) close to 1, to the data for all system sizes (L > 6 ). 
These allow us to conclude that g c = 4.83704(6), where 
the numbers in parenthesis denote the expected errors 
(one standard deviation) in the preceding digit, i.e. the 
relative error on g c is approximately one part in 10 5 . If 
we allow independent parameters gf} and g£ in the fits to 
the two data sets, both estimates of the critical point are 
statistically consistent with this g c , albeit with somewhat 
larger error bars. 

We note here that the values we find for the subleading 
exponent, u = —0.31(5) for the R 2 data and lo = 0.82(5) 
for the p s L 2 data, he far from a common asymptotic 
value. Thus indeed u should be considered as an effective 
exponent accounting for crossover effects in system size, 


The critical O(N) </> 4 theory, by which is meant the 
theory at the upper critical dimension and at the criti¬ 
cal point, obeys many fundamental and universal prop¬ 
erties, some of which depend on N while others are 
IV-independent. In Ref. |35| it was shown that the ze¬ 
ros of the partition function (Lee-Yang zeros) [35], and 
hence the thermodynamic functions, obey a finite-size 
scaling theory, which was derived by renormalization- 
group methods. These perturbative arguments demon¬ 
strate that there are multiplicative logarithmic correc¬ 
tions in the system-size dependence of derivable thermo¬ 
dynamic functions, which are closely finked to those of 
the Lee-Yang zeros and, furthermore, are independent of 
N for odd N. This leads to the key practical observa¬ 
tion that size-dependent logarithmic corrections in phys¬ 
ical observables such as the magnetic susceptibility and 
the specific heat at the critical point follow a universal, 
IV-independent form when N = 3. Here we provide a 
non-perturbative calculation of the magnetic susceptibil- 
ity, x(Qaf ,L) [Eq. ([ 8 ])], for systems of finite L at the 
QCP, g c , of the (3 + l)-dimensional 0(3) transition to 
test the predicted logarithmic corrections. 

The universal form of the magnetic susceptibility at 
the critical point in a finite-size system is given by 13? 


%(Qaf ,L) = aL 2 [lnL ] 1/2 


1 + 6 


ln(lnL) 

InL 


(19) 


with non-universal but L-independent parameters a and 
6 . We used this expression with a fixed value g = 4.837, 
which is within the standard deviation of the g c value 
found in Sec. m to investigate the logarithmic cor¬ 
rections to the L-dependence of x(Qafj^)- We cal¬ 
culate the susceptibility at the ordering wave vector, 
Qaf = (7T, 7T, 7r, 7r) , for systems of all even sizes from 
L = 6 to 40. As in Sec. m the scaling predictions un¬ 
der test are valid for a 4D theory and again we use the 
critical contour T -1 = kL with k = 1. 

Figure [3|a) shows our results for the critical mag¬ 
netic susceptibility normalized by L 2 . In the absence 
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FIG. 3. (Color online) (a) QMC data for x(Qaf,£)/£ 2 ob¬ 
tained at g = 4.837 for all even system sizes from L = 6 to 40. 
Solid lines are fits to a(lnL) 1 ^ 2 (green) and to Eq. ( |l9| ) (blue). 
We apply the square-root fit only for system sizes L > 30 and 
the optimal value of the fitting parameter is a = 0.274. The 
two-parameter fit is made to the data for all system sizes 
L > 14 and yields optimal parameters a = 0.522, b = —1.317. 
(b) Reduced y 2 values obtained by fitting \'(Qaf ,L)/L 2 , for 
14 < L < 40, to the form ( |19| ) , but with the exponent 1/2 
replaced by a parameter fj. The optimal y 2 value is obtained 
at fj ~ 0.5, consistent with the prediction of Ref. [35]. 


of logarithmic corrections, x(Qaf ,L)/L 2 would be con¬ 
stant and the curve would be a flat line. Instead we 
observe that a reasonable account of the data for our 
larger system sizes (L > 30) requires a fit of the form 
X (Qaf ,L)/L 2 = a(lni) 1 / 2 , as anticipated in Ref. [35] . 
However, for a more quantitative fit over the full size 
range available, we find [Fig. |3][a)] that it is necessary to 
include the predicted subleading logarithmic correction 
in Eq. (19). 

To examine the sensitivity of these results to the ex¬ 
ponent 1/2 of the multiplicative logarithm in Eq. (19), 
we replace this predicted exponent by a variable fj. We 
determine this exponent by calculating the goodness of 
fit Xr as a function of fj. As Fig. [3](b) makes clear, the 
best fits are indeed obtained close to fj = 0.5, in complete 
consistency with Eq. (19). 



FIG. 4. (Color online) Staggered magnetization, defined for 
each system size as m s (L) = [3(m2 Z (T))[ 1 ^ 2 , shown as a func¬ 
tion of 1 /L 2 for a range of coupling ratios g near g c . Poly¬ 
nomial fits of cubic order were used to extrapolate m s (L) to 
the thermodynamic limit; the temperature in all cases was 
T = 1/L. Error bars on all points are smaller than the sym¬ 
bol sizes. 


From the fact that our exact numerical data confirm 
not only the leading but also the subleading corrections 
to scaling, we conclude that obvious logarithmic cor¬ 
rections can be observed in the size-dependence of the 
thermodynamic functions at the QCP. This result also 
demonstrates that our determination of g c is sufficiently 
precise to study logarithmic corrections without signifi¬ 
cant distortions arising from uncertainties in its value. 


V. SUBLATTICE MAGNETIZATION 

Physical condensed-matter systems at continuous 
QPTs are generally in the thermodynamic limit, and 
size-scaling measurements of the type easily performed in 
QMC simulations (Sec. IV) are not a realistic experimen¬ 
tal option. However, as discussed in Sec. I, multiplica¬ 
tive logarithmic corrections are expected in a range of 
physical quantities close to the QCP. The primary phys¬ 
ical observables in the quantum antiferromagnet are the 
zero-temperature staggered magnetization, m s (g c — g ), 
and the Neel temperature, Tpj{g c — g). Calculating 
these quantities in the thermodynamic limit is signifi¬ 
cantly more challenging than studies of size-dependence, 
as careful extrapolations of finite-size data are required. 
Here and in Sec. VI we describe and then implement ap¬ 
propriate measures for extrapolating to infinite system 
size and, for m s , to zero temperature, thereby revealing 
the logarithmic corrections to both m s and T (y. 

We compute the staggered magnetization according to 
Eq. 0 for coupling ratios as close to g c ~ 4.837 as 
g = 4.834. For a given value of g , we calculate the 
squared quantity (m^ z (g, L )) over a range of system sizes. 




















As shown in Fig. [4j the staggered magnetization clearly 
decreases with increasing L and converges towards a fixed 
limit, suggesting a controlled extrapolation even very 
close to the QCP. Definitive extrapolation to the ther¬ 
modynamic limit in this regime is a complex issue, and 
a discussion of several technical points is in order before 
analyzing our results. 


A. Extrapolation scheme 

First, most of the simulations in this section are per¬ 
formed at a temperature T = 1 /L, such that the ex¬ 
trapolation L —> oo includes both system size and tem¬ 
perature. Because the system is ordered for sufficiently 
low T and the order parameter converges quickly to a 
non-zero value below the ordering temperature, one may 
use the form T = aL~ b with arbitrary prefactor a > 0 
and exponent 6 > 0 to study the T —>■ 0 magnetiza¬ 
tion as a function of L. These conditions for obtaining 
a (T — > 0,P — > oo) extrapolation of the order parameter 
for g < g c contrast with the need to follow the contour 
T = aL~ 3 (a = 1/fc; 6=1 when z = 1) in Secs. Ill 
and IV for studies of the QCP. Although large values 
of a and 6 should improve the convergence, in practice 
one must consider the balance between computation time 
and convergence rate, and the choice a = 1 , 6 = 1 works 
well in most cases. However, for coupling ratios very 
close to the QCP, the temperature may be a significant 
fraction of Tjv and thus m s (g : L) could be far from its 
zero-temperature value. We have therefore performed 
additional simulations at T = 1/(2L) to verify that the 
extrapolation does remain well controlled and fully rep¬ 
resentative of the thermodynamic limit in temperature 
as well as in system size. 

Second, in contrast to Sec. a where we used the non¬ 
trivial power-law scaling forms (161 known to be appro¬ 


and ( 21 ) will not be exactly the same, but for reliable fits 


priate for extrapolating the location of a critical point, 
the ground-state order parameter inside the Neel phase 
can be extrapolated by using simple polynomial fits. To 
obtain m s , one may extrapolate the squared quantity 
and then take its root afterwards or take the square root 
for each system size before extrapolating (the procedure 
followed in Fig. [4]); the corresponding polynomials are 

™?sz(g , L ) = a (g) + Kg) L ~ 2 + c(g)L~ 3 +..., ( 20 ) 
Vm 2 S z{g,L ) = a'{g) + b'(g)L~ 2 + c'(g)L ' 3 + ...(21) 

Here the leading L-dependence in the extrapolation of 
a non-vanishing order parameter, at fixed g inside the 
ordered phase, is known [32] to be L 2 ~ d ~ z due to the 
dimension-dependent power-law decay of the transverse 
correlation function. Here d + z = d + 1 = D = D c = A, 
and the resulting leading 1 /L 2 dependence is shown 
clearly in Fig. [3] Because these procedures use a poly¬ 
nomial of finite order to approximate physical behavior 
containing, in principle, an infinite number of corrections, 
the extrapolated value of m s obtained with the forms ( 20 ) 


they should agree within statistical errors. 

We stress that no logarithmic corrections are expected 
in this case, meaning that on grounds of principle they 
should not be present in the asymptotic large-L correc¬ 
tions to a non-zero-valued order parameter. This non- 
critical behavior contrasts with the case of the shift in 
the critical point discussed in Sec. m where logarithmic 
corrections should in principle be present, although we 
concluded that their effects are not detectable in prac¬ 
tice. 

Finally, however, non-trivial corrections may still be 
expected in the L dependence of ( m 2 z (g,L )) for g close 
to the QCP, where the order parameter is small, in the 
form of crossover behavior from near-critical at small sys¬ 
tem sizes to asymptotic ordered-state scaling at large L. 
Quite generally, no analytic functional forms are available 
for describing such crossovers, and great care is required 
to ensure that the asymptotic region, where Eqs. (20) 
and ( 21 ) are valid, has been reached. As g -A g c , succes¬ 


sively larger systems are required for this, and here we 
find that reliable extrapolations are no longer possible 
beyond g = 4.834 because of the limits on system size 
set by the available computer resources. 


In fits to the forms (20) and (21), it is necessary to 


select the order P of the polynomial and the range of 
system sizes to include. The size of the error bars on the 
QMC data points has a significant influence here, because 
deviations from the leading L -2 correction are easier to 
detect with smaller error bars. We have performed a 
systematic study using fits of orders P = 3 to 6 , including 
different system-size ranges. We characterize the quality 
of the fits using the standard reduced % 2 measure, and 
for a “good” fit we require that the optimal value must 
fall within three standard deviations of its mean, i.e. we 
demand that 


Xr - 1 = 


X 


- 1 < 3, 


n L - n p 


n L - n p 


( 22 ) 


where is the number of data points (system sizes) and 
n p = P + 1 is the number of fitting parameters. For a 
given P and largest system size P, we use all available 
system sizes down to a smallest size L m - m for which the 
above condition is still satisfied. We then study the be¬ 
havior as a function of L for different P, and compare 
the values of m s obtained from extrapolations based on 
Eqs. (20) and (21). To estimate the error bars on the ex¬ 
trapolated m s (g), we performed additional polynomial 
fits with Gaussian noise (whose standard deviation is 
equal to the corresponding QMC error bars) added to 
the finite-size data. The standard deviation of the distri¬ 
bution of extrapolated m s (g) values defines the statistical 
error. 

Figure [5] shows results for several values of g approach¬ 
ing the QCP, with cases where the fits are relatively 
straightforward (further from g c ) shown in panels (a) and 
(b) and more challenging cases (closer to g c ) shown in 
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FIG. 5. (Color online) Extrapolated values of the sublattice 
magnetization as a function of the largest system size included 
in the fit, for values g near the QCP ( g c ~ 4.837). Panels 
(a) and (b) include values of g for which the extrapolations 
are stable for relatively small L; panels (c) and (d) show g 
values very close to g c , where the extrapolations require large 
sizes to stabilize. Results for different orders P of the fitting 
polynomial are compared. The smallest system size included 
was determined using the \r criterion of Eq. (22). Panels 


(a,c) and (b,d) show respectively the results of extrapolations 
of (m 2 z (g, L)) and {m 2 z (g, L)) 1 ' 2 , with the square root taken 
after the extrapolation in the former case. 


panels (c) and (d). The upper panel in each group corre¬ 
sponds to the square root being taken after the extrapola¬ 
tion [Eq. (20)], while the lower corresponds to fitting the 
square root for each system size [Eq. (|2l])]. In panels (a) 
and (b), the extrapolated values are observed to be very 
stable with respect to the range of system sizes and the 
order of the polynomial, whereas panels (c) and (d) man- 



FIG. 6. (Color online) Extrapolated staggered magnetization, 
m s at T = 0, as a function of the distance from the QCP, 
using the value g c = 4.83704 determined in Sec. |III[ Error 
bars on the calculated data points are similar to or smaller 
than the symbol size. The closest point to g c is g = 4.834. 
Lines show both the best fit by a pure square-root function 
[Eq. ( |23| ), green] and including the logarithmic correction fac¬ 
tor predicted in Ref. [35] [Eq. (241, blue]. The fitting pa¬ 
rameters of the logarithmic correction curve are a = 0.266(2) 
and b = 4.8(3). The yellow shading represents the approxi¬ 
mate extent of the QC regime and is determined by including 
all data points described adequately (within a deviation of 
approximately 4%, see text) by the functional form of the 
logarithmic correction curve. The inset shows rn s {g) and the 
QC regime on linear axes. 


ifest some of the crossover behavior expected close to g c , 
showing considerable variation as the maximum system 
size is increased. There are also significant differences 
between the two fitting procedures, until the largest sys¬ 
tem sizes where the extrapolations stabilize; we take the 
fact that the two types of fits give consistent results for 
these largest systems at all values of g as an indication 
that the extrapolations are reliable. We have not been 
able to achieve good convergence based on system sizes 
up to L = 48 for g values closer to g c than those shown in 
Figs. [5][c) and (d). In these most challenging cases, our 
results show that it is better to use the fitting form of 
Eq. extrapolating the square root of the staggered 
magnetization for each system. Regarding the quality of 
the fits obtained by varying the polynomial order, P, in 
Fig.0 we find that extra terms in the fit scarcely justify 
the additional degrees of freedom lost in the determina¬ 
tion of Xr- All of the results presented below were ob¬ 
tained by extrapolating yjm 2 with polynomials of order 
P = 4. 


B. Thermodynamic Limit 

With all of the above considerations, we are able to 
obtain reliable and high-precision extrapolations of the 
staggered magnetization in the thermodynamic limit for 
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values of g as close to the QCP as |g — g c \ ~ 0.003. In 
Fig. [ 6 ] we show all of our data for m s (\g — g c \) on log¬ 
arithmic axes. If these data satisfied mean-field scaling 
alone, with no discernible logarithmic corrections, one 
would expect a curve of the form 

ms(g) = a\g - g c \ 1/2 , (23) 

but this (green line in Fig. [ 6 ]) is manifestly unable to de¬ 
scribe the data. For the zero-temperature order param¬ 
eter, perturbative renormalization-group considerations 
applied to the O(N) cf) 4 field theory at the upper critical 
dimension predict the form 

m s (g) = a\g - g c \' 3 \ln(\g - g c \/b)f , (24) 


where /3 = 1/2 is the mean-field exponent and the expo¬ 
nent of the multiplicative logarithmic correction is given 
by j3 = 3/(N + 8) [33- A fit to this form, using /3 = 3/11 
for N = 3 (blue curve in Fig. [b|, yields excellent agree¬ 
ment with the data all the way to our smallest values of 
\g — g c \\ the fitting parameters are a = 0.266 ± 0.002 and 
b = 4.8 ± 0.3. We note that the fit is very insensitive to 
the precise value of b , and for further analysis we fix this 
to b = g c . 

To test the predicted exponent $ = 3/11 in Eq. (241, 
we treat it as a free parameter and fit our data using 
different numbers of g values, including all points closest 
to g c and studying the behavior as points further away 
from g c are added one by one. Figure [7] shows y 2 and /3 
as functions of the number of data points fitted. With 
the exception of cases including the two points furthest 
away from g c , all the fits appear reasonable, with % 2 < 2 . 
However, by the properties of the y 2 distribution, a fit 
should be considered statistically acceptable only if a cri¬ 
terion analogous to Eq. ( [22] ) is satisfied, i.e. the largest 
number of data points for which y 2 — 1 remains less than 
three times its standard deviation (3cr) marks the bound¬ 
ary between good and poor fits. At this point we obtain 
j3 = 0.268 ± 0.008, which lies well within one standard 
deviation of the predicted value 3/11 ~ 0.2727. If more 
points are excluded, the fitted exponent evolves slowly 
[Fig. [7|b)] while remaining statistically well compatible 
with the predicted value. Because the fitting error in¬ 
creases, less weight should be placed on results including 
less data, and taking an error-weighted average over all 
the points below the cut-offline, N g = 23, in Fig. [TJyields 
j3 = 0.279±0.011. We take this as complete confirmation 
of the predicted value. 

As important as finding clear logarithmic corrections 
to scaling is that we have demonstrated their presence 
over a significant region around the QCP; indeed, most 
of the points we have computed are well described by 
Eq. (24). Including the multiplicative logarithmic correc¬ 
tion converts an inadequate description of the data into 
an excellent one (Fig. [6]) as far inside the Neel phase as 
\g — g c \/9c ~ 0 . 2 , where the order parameter is already at 
60% of its maximum possible value (m s = 1 / 2 , at which 
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FIG. 7. (Color online) Statistical an alys is of the exponent 
of the logarithmic correction in Eq. (241. (a) Reduced y 2 

value of the fit, normalized to the standard deviation, and 
(b) optimal value of the exponent shown as functions of the 
number of data points (g-values) used, beginning from the 
point closest to g c in Fig. [6] The vertical dashed line indicates 
the numbe r of points, N g = 23, for which a 3<r criterion for 
Xr [cf. Eq. (22 i] is satisfied, as indicated by the horizontal line 
in panel (a). In panel (b), the error bars were computed by 
repeating the fits multiple times with Gaussian noise added to 
the m s data points. The horizontal line marks the predicted 
value /3 = 3/11. 


point no quantum fluctuation effects remain). This im¬ 
provement is clearer still in the inset of Fig. [ 6 j which 
shows the results on linear axes. Under the assumption 
that data points at large \g — g c \ no longer fall on the fit¬ 
ted curve because they lie outside the region controlled 
by the QCP, we can determine the size of the critical 
region based on a threshold maximum deviation of the 
data from the curve. Although the choice of threshold 
value is somewhat arbitrary, the |g — g c \/gc < 0.2 region 
indicated by the yellow shading in Fig.[4]reflects a thresh¬ 
old of approximately 4%, which lies well above achievable 
experimental uncertainties. We comment in Sec. |VII| on 
the utility of our results for the case of TlCuCU. 
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VI. NEEL TEMPERATURE 

We turn next to the scaling form of the Neel temper¬ 
ature, Tjv(fl'), near the QCP. Unlike the T = 0 order pa¬ 
rameter, as far as we are aware there is no prediction from 
perturbative field-theoretical calculations including loga¬ 
rithmic corrections for the scaling form of finite-T critical 
points at the upper critical dimension. Close to a QCP, 
the general power-law form without logarithmic correc¬ 
tions is discussed in Ref. 3], but one may also expect a 
multiplicative logarithmic term as in the other quanti¬ 
ties we have discussed. We first derive the exponent of 
the logarithm for the 0(3) transition in 3+1 dimensions, 
based on the known scaling properties of related quan¬ 
tities. We then present our QMC calculations of T/v for 
the Heisenberg model on the double cubic lattice and test 
our prediction. 


A. Scaling hypothesis 


In a path-integral construction in imaginary time, the 
size of the system in the time dimension is proportional 
to the inverse temperature, /3 = 1/T. This can be con¬ 
sidered as a length, T t , on a parallel with the spatial 
lengths, T, of a d + 1-dimensional system. If the spa¬ 
tial lengths are taken to infinity in all directions, what 
remains is a single finite length, L T , for the effective sys¬ 
tem, and finite-size scaling in this length corresponds to 
finite-T scaling in the original quantum system 2Til . 

Without logarithmic corrections, by analogy with the 
finite-T shift of the critical point discussed in Sec. ME 
the same type of shift as in Eq. (15) can be expected 
because z = 1. Thus 


9c(T) - g c { 0) ~ L- 1 '”, 


(25) 


as a consequence of the finite temporal size, and the scal¬ 
ing behavior is Tjv ~ (g c — g ) v , as discussed in detail in 
Ref. p). In the case of spatial finite-size scaling, with all 
lengths finite, the shifted critical point (sometimes called 
the pseudo-critical point) is not a singular point, but the 
singularity develops as L —> oo. By contrast, in the finite- 
T case in d = 3 spatial dimensions, the shifted point is a 
true (classical) phase transition, although from a scaling 
perspective this difference is not relevant. 

In order to discuss logarithmic corrections, it is useful 
to first express T/v using a macroscopic, zero-temperature 
energy scale of the system that vanishes as g —> g c [3]. 
For the spin system considered here, the only such energy 
scale is the spin stiffness, p s . According to Ref. I3Z1. the 
scaling form of this quantity in the ordered phase when 
z = 1 is 


p s ~(5 c-ffr (d - 1} . (26) 

Consistency with the result T/v ~ ( g c — g) v then gives 
the scaling of the critical temperature for d = 3, 


(27) 


where the mismatch in units is compensated by a power 
of the non-singular spin-wave velocity [3] (Sec. VIC), 
which can be neglected here. Our basic hypothesis is 
that this proportionality, which is the singular part of 
a relationship based on matching scaling dimensions, ap¬ 
plies in all respects at the upper spatial critical dimension 
(d = 3 for 2 = 1), such that logarithmic corrections to T n 
arise solely due to the logarithmic corrections intrinsic to 
Ps- 

Fisher et al. [221 have shown that the critical spin stiff¬ 
ness can be expressed as p s ~ £ 2 /, where £ is the corre¬ 
lation length and / is the free-energy density. The loga¬ 
rithmic corrections to both £ and /, presented by Kenna 
in Ref. [25], are 

\g~ 9 c\~ 1 ' hf(\g -g c \), (28) 


with v = 5/22 for the relevant universality class, and 

/ ~ Iff — 3c| 4i/ ln “(|<7 — <7c|). (29) 

with a = 1/11. The logarithmic correction to p s is there¬ 
fore given by 


p s ~|s-Sc| 2 i 5 ln 2 i>+a (|s-Sc|), 


(30) 


and by combining these results with Eq. (27) we obtain 


Tn ~ |s-Scrin T (|s-Sc|), 


(31) 


where f = v + d/2. From the values of v and a given 
above [25], we obtain the prediction f = 3/11, which 
is remarkable in that the exponent in the logarithmic 
correction to TV should be the same as the one in the 
sublattice magnetization, t = j3 (24). 


Because the zero-temperature order parameter is a 
consequence purely of quantum fluctuations, whereas the 
classical ordering temperature is a consequence primar¬ 
ily of thermal fluctuations, there is a priori no reason 
to expect that the two should have the same form. Ex¬ 
act numerical calculations are therefore uniquely posi¬ 
tioned to provide qualitatively new information in this 
case. We note that this equality applies to the phase 
transitions of O(N) models for all values of N; because 
v = (V + 2)/[2(7V + 8)] and a = (4- N)/(N + 8) gg, we 
obtain f = 3 /(N + 8), the same value as the exponent /3 
in Eq. (24). Thus we predict that, in the neighborhood 
of g c , Tiv(g) will be proportional to m s (g,T = 0), with 
no multiplicative logarithmic factors, for all values of N; 
this result was reported for the N = 3 case in a previous 
QMC study [5D], which we now extend sufficiently close 
to the QCP to observe the cancellation of logarithmic 
terms. 


B. QMC calculations 

Calculating Tfj(g) within our QMC simulations is sim¬ 
ilar to obtaining g c in Sec. |III[ but with some impor¬ 
tant differences of detail. The calculations of Sec. IIIII 


rj~ i 2 

T/v ~ Ps, 
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were performed for a genuinely 4D system, with the 
imaginary-time axis treated on the same footing as the 
spatial dimensions. At finite temperatures, this symme¬ 
try is broken and the system is 3D with a separate tem¬ 
perature variable, which determines the finite thickness 
of the time dimension even when L —> oo. Both the 
Binder ratio [Eq. ([£])] and the spin stiffness [Eq. ©I are 
size-independent quantities at the thermal phase transi¬ 
tion and hence remain valuable indicators, although the 
appropriately scaled spin stiffness for the 3D transition 
is now p s L (instead of p s L 2 , used for analyzing the 4D 
T = 0 transition in Sec. Ill I | i2U] . 


For each value of the coupling ratio g within the Neel 
phase (g < g c ), we compute R 2 and p s L for a range 
of system sizes and perform finite-size-scaling extrapo¬ 
lations to deduce the Neel temperature, Tiy(g), in the 
thermodynamic limit. Similar to Sec. |III[ we first ob¬ 
tain the crossings of the R. 2 {T) and p s L(T) data for dif¬ 
ferent system sizes using polynomial fits, as shown in 
Figs. [8[a) and[8[b) for g = 4.71. The crossing points of 
both quantities for each successive pair of system sizes, 
T N (g,L ) and T N (g,L + 2), are used to extrapolate to¬ 
wards the value T/v(g,L —>• 00 ) from above and below, 
using power-law forms analogous to Eq. (16). We note 
that, as in the analysis leading to g c (Sec. Ill B), the 
data points obtained for 7? 2 (A,T) and p s L(L,T) using 
systems of all sizes (L > 10) fall within a 3<r criterion 
analogous to Eq. ( [22] ) for this type of fit. The extrapola¬ 
tion of T]y(g, L —» 00 ) = Tjv(<?) = 0.5363(13) for g = 4.71 
is shown in Fig. |8]c). 


We comment here that our determination of Tjy(g) for 
g close to g c is rather less precise than our determination 
of m s {g). The fundamental difference in character of 
the two quantities, and hence of their calculation, causes 
the estimators for Tjy(g) [the approximate crossings in 
Figs. |8]a) and[8[b)] to have larger error bars and finite- 
size effects. Further, the error bars of the crossing points 
grow rapidly as g —> g c , while the decrease in T)v leads to 
longer simulation times (because the space-time volume 
is proportional to L 3 /T). After detailed error control, 
the closest reliable data point to the QCP is g = 4.831, 
for which \g — g c \ is twice as large as for the closest m s (g) 
point (Sec. [V]. We have nevertheless obtained 19 reli¬ 
able data points, within the QC regime determined from 
m s (g) (Fig. [6]) and down to unprecedentedly low temper¬ 
atures, which are fully sufficient to test for evidence of 
logarithmic corrections to T/v(g). 

The Neel temperature has units of energy and clearly 
depends on the overall energy scale of the system. Ideally, 
it should be normalized by an intrinsic energy scale of the 
system to give a dimensionless quantity. In Ref. [HI it 
was shown that T^{\g — g c |) normalized to the micro¬ 
scopic energy scale J s , given by the sum of all couplings 
of a spin to its neighbors, yields a remarkably system- 
independent result for Heisenberg antiferromagnets with 
three different dimerization patterns; for the double cu¬ 
bic lattice, J s = J(6 + g). Other authors EH [22] have 
suggested that the appropriate normalization is given by 




1/L 


FIG. 8. (Color online) Procedures used to extract the Neel 
temperature, Tjv(<?), illustrated for the case g = 4.71. (a) 
Binder ratio R 2 as a function of T for system sizes L = 10, 
12, ..., 30. (b) Scaled spin stiffness p a L as a function of T for 
the same values of L. Error bars are smaller than the symbol 
sizes. Crossings of these lines are extracted using polynomial 
fits and the results are used to obtain finite-size estimates for 
quantities T$(L) and T^(L). (c) Fits to data for systems of 
all sizes (L > 6) of the two size-dependent crossing estimators 
using functions of the form Tjv(L) = Tjv(oo) + a/L 1/,1 ' + “. 
Enforcing the same constant Tjv(oo) = Tjv for the R 2 and p s L 
crossings gives Tn = 0.5363(13) with irrelevant exponents 
to « 0.8(2) for R 2 and 1.1(3) for p s L (1/v « 1.42 for the 
relevant 3D universality class). 


a macroscopic quantity, the spatially averaged spin-wave 
velocity, ^/c x c y c z (which, it should be noted, does not 
have units of energy and requires an unknown dimension¬ 
ful constant). We begin by taking the former approach 
and return below to address the latter. 

The relationship between Tn/J s and \g — g c \ is pre¬ 
sented in Fig. [9] Once again we show a mean-field scal¬ 
ing line for comparison and once again it cannot provide 
an adequate fit, suggesting that logarithmic corrections 
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FIG. 9. (Color online) Normalized Neel temperature, T^/J s , 
as a function of the distance from criticality, |g — g c \. The 
closest point to g c « 4.837 is g = 4.831. Lines show both the 
best fit by a pure square-root function (green) and using a log¬ 
arithmic correction factor with exponent f = 3/11 [Eq. (311, 
blue]. The yellow shading represents the QC regime and is 
determined from Fig. [5] The inset shows Tjv(g) on linear axes. 


are indeed present. However, a fit to our predicted form, 
given by Eq. (31) with f = 3/11, describes the data very 
well, even at the limits of the region classified as QC 
based on the m s (g c — g) fit in Fig. [5] 

For a fully quantitative test of the exponent we predict 
for the multiplicative logarithmic term in Eq. (31), we 
substitute a free exponent f for the fixed value 3/11 and 
optimize it using fits with different windows of g- values. 
This analysis is precisely analogous to that performed for 
m s in Fig. |T] The behavior of Xr and °f the optimized ex¬ 
ponent, with error bars again estimated using the method 
of numerical Gaussian noise propagation, is presented in 
Fig. |TU} By taking the inverse-variance-weighted aver¬ 
age over all results for which Xr is acceptable, we obtain 
f = 0.275(2), in excellent agreement with the prediction 
f = 3/11 ~ 0.2727. We conclude that the multiplicative 
logarithmic correction to T^(g) is, to within our error 
bars and in agreement with a straightforward scaling ar¬ 
gument based on the spin stiffness (Sec. VIA), identical 
to the m s {g) correction in Eq. ([24} . 


C. Spin-wave velocity 


The spin-wave velocity, c, is uniform in the primary 
axial directions on the double cubic lattice. As discussed 
in Sec. |nj it can be calculated most straightforwardly and 
most accurately in the SSE framework from the spatial 
and temporal winding-number fluctuations in Eq. (13) 
to define the space-time-isotropic criterion of Eq. (14), 


which contains the velocity PH [22} HU [46]. This tech 
nique remains well-defined throughout the critical regime 
and is expected not to be affected by logarithmic cor- 



Number of Points 


FIG. 10. (Color online) Statistical a naly sis of the exponent 
of the logarithmic correction in Eq. (311, performed by re¬ 


placing the predicted value 3/11 with a fitting parameter r. 
(a) Reduced y 2 value of the fit, normalized to the standard 
deviation, and (b) optimal value of the exponent, both shown 
as functions of the number of data points (g-values) used, 
beginning from the point closest to g c in Fig. [9] The ver¬ 
tical dashed line indicates the number of points, N g = 16, 
included in the fit below which y 2 satisfies a 3<r criterion anal¬ 
ogous to Eq. (221, as indicated by the horizontal line in panel 
(a). In panel (b), the error bars were computed by repeating 
the fits multiple times with Gaussian noise added to the TV 
data points. The horizontal line marks the predicted value 
t = 3/11. 


rections; it was shown in Ref. [16] . which we follow for 
technical details, that the winding-number approach pro¬ 
duces the correct result for c in the Heisenberg chain, a 
system known to have strong logarithmic corrections to 
scaling. 

In Fig. 11 a) we show the results for c(g,L) of cal¬ 
culations on finite systems of even sizes up to L = 26, 
which we extrapolate to the thermodynamic limit using 
the relation 


c(g, L ) = c(g) + a{g)/L 2 + b(g)/L 3 . (32) 

This form is found empirically j!6j to provide a very good 
reproduction of the data and numerical errors due to 
finite-size effects in the critical regime are clearly small. 
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l/L 



FIG. 11. (Color online) Calculation of the spin-wave velocity, 
(a) Velocities c(L) obtained for systems of sizes from L = 6 
to 26 at different values of the coupling g and extrapolated 
using Eq. ( |32[ ). (b) Extrapolated velocities c(g) as a function 
of the microscopic energy scale 6 + g. Error bars are mostly 
hidden inside the symbols. The solid line is a linear fit and 
the inset magnifies the region close to g c . 

Figure [TTJb) shows the results for the extrapolated spin- 
wave velocities c(g) of the infinite system, by comparing 
the macroscopic scale c 3 ^ 2 with the microscopic quan¬ 
tity 6 + g discussed above. The almost perfect linearity 
demonstrates that the two effective energy scales are very 
closely related, which can be expected from the fact that 
the velocity of a spin excitation depends directly on the 
net interaction of a single spin, and both are perfectly 
valid choices for the normalization of T/v- A graph com¬ 
pletely analogous to Fig. [9j showing a logarithmic cor¬ 
rection with the same exponent f = 3/11, is obtained if 
c 3//2 is used to normalize T/v(g). 


D. Relation between Tn and rn H 

In experiment it is often difficult to relate an external 
control parameter to the microscopic coupling constants 
of a model Hamiltonian. In a quantum antiferromag- 
net, some aspects of this problem can be circumvented 
by studying the relationship between Tn and m s (T = 0 ) 


directly, without reference to the control parameter, g. 
A universal relationship between these macroscopic and 
measurable quantites would be of considerable experi¬ 
mental utility in characterizing the nature of critical phe¬ 
nomena without recourse to detailed microscopic knowl¬ 
edge of the system parameters (such as the pressure de¬ 
pendence of the exchange couplings in TICUCI 3 ). Al¬ 
though an experimental test [ 8 ] of the linear relationship 
between Tn and m s [20j indicated satisfactory agreement 
close to the QCP, the issue of how best to normalize Tn 
was not addressed. We use our systematic data spanning 
the entire QC regime to test the limits of linear propor¬ 
tionality and discuss the normalization of Tn- 

In Ref. 1201 , where a universal linear relation was found 
in three different models, the authors articulate a mean- 
field argument based on semiclassical considerations for 
a direct proportionality of Tn to the effective spin order 
gauged by m s at T = 0. In Ref. [2T], these arguments 
were elucidated in a field-theory context, where it was 
stated that logarithmic corrections should be negligible 
for linear proportionality to emerge. In fact these argu¬ 
ments can be reduced to the statement that it should be 
possible to treat quantum and thermal fluctuations in¬ 
dependently, with no mutual interference of their effects 
[2Uj . If one considers that mean-field exponents are valid 
in high-dimensional systems ( D > D c ) because ther¬ 
mal fluctuations become independent of quantum fluc¬ 
tuations when the phase space is sufficiently large, then 
it appears that weak logarithmic corrections could enter 
the relationship of Tn to m s at D = D c . This possibility, 
also motivated by the (then) unknown form of the log¬ 
arithmic corrections to Tn, was investigated directly by 
QMC simulations for the cubic lattice [22] , but the results 
were not conclusive (claims concerning the observation of 
logarithmic corrections are not justified by the available 
data range). Here we have presented scaling arguments 
(Sec. VIA) and numerical data (Sec. VIB) demonstrat¬ 
ing that the logarithmic corrections to m s and Tn have 
precisely the same form, setting their linear relationship 
in this class of system on a far firmer foundation. 

Our data from Secs. [V] and [Vi] can be used to probe 
the T)v(to s ) relation in detail and confirm that linearity 
extends much closer to the QCP than previous studies 
could show. Because our results are for a single type of 
dimerized model, we are not able to address the question 
of a universal prefactor [2U]. However, we are able to 
make a definite statement regarding logarithmic correc¬ 
tions in the relationship between Tn and m s . Figure [l2| 
shows our data, taken from Figs. [ 6 ] and [9] in the form 
X)v(to s ), with the implicit control parameter g effectively 
eliminated. In panel (a), Tn is simply normalized by the 
energy scale, the interdimer coupling J = 1; in panel (b), 
we have normalized Tn by the composite scale J s = 6 + g 
(where g can be considered a function of m s ), as in Fig. [9j 
in panel (c), we have normalized Tn using the correctly 
scaled spin-wave velocity, c 3 / 2 . The shaded regions again 
signify our definition of the critical region, based on the 
strict critical scaling of rn s (g) in Fig.[5j 
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m s 


FIG. 12. (Color online) Relationship between Tjv and iris, 
using different normalizations of T)v; in (a) by the intra-cube 
coupling J = 1, in (b) by the sum of couplings J a = 6 +g, and 
in (c) by the spin-wave velocity in the form sj 12/hc 3 ^ 2 . The 
lines are linear (proportionality) fits to the small-m s points 
and the yellow shaded areas denote the QC regime determined 
from m s (g) in Fig. [5] 


As shown in Ref. [20], we find in Fig. 1 1 2[b) that the 
linearity of m s and Tn/J s extends well beyond the QC 
regime; although our data were not selected to focus on 
this region, we do find complete agreement with previ¬ 
ous calculations where the data overlap. Here we demon¬ 
strate that the essentially perfect linearity also extends 
much closer to m s = 0 = T'n (g g c ), and in particular 
that it remains valid throughout a regime with explicit 
logarithmic corrections in the individual quantities. Al¬ 
though we cannot show that the linear relationship ex¬ 
tends all the way to the QCP, our data certainly suggest 
that this is the case, i.e. that it can also be considered 
a universal property of the QC regime. To the extent 
that linearity of m s and T'n is a consequence of the de¬ 
coupling of thermal and quantum fluctuations, this inde¬ 
pendence appears to extend from strongly ordered sys¬ 
tems, where no logarithmic corrections are expected, to 
the most strongly fluctuating QC systems. The linear re¬ 
lationship we have demonstrated verifies in full our scal¬ 


ing prediction that the logarithmic corrections to T\ r (g) 
have the same exponent, f = $ = 3/11, as m s {g). We 
also note that the normalization of T)v has no effect (be¬ 
yond the prefactor) on the linear relationship in the QC 
regime, but that deviations from linearity clearly differ 
at higher m s . 

We close by considering in more detail the case where 


T)v is normalized by c 3 / 2 [Fig. 12 [c)], with a view to 


making quantitative comparisons with field-theory pre¬ 
dictions [21] . The explicit relationship is 


T n = 7 c' 


.3/2 


(33) 


where 7 = (ip) /m s is the dimensionful prefactor relating 
the expectation value of the un-normalized field /> in the 
action to the order parameter m s of the lattice model. 
In the dimensionless units of our work (J = 1 ,K = 1, 
...), we obtain 7 = 0.6998 ± 0.0016, thereby providing a 
bridge between the quantum field theory and the micro¬ 
scopic lattice Hamiltonian. We suggest that this calcu¬ 
lation should be repeated for other dimerized geometries 
to test the universality of 7 . It is worth repeating in 
this context the advantages of the double cubic lattice in 
making the spin-wave velocity equal in all three primary 
axial directions. The winding-number method used to 


extract c in Sec. VIC can be generalized to anisotropic 
systems E3, but incurs the significant complication of 
altering the aspect ratio of the spatial lattice. Different 
techniques for computing the velocities, such as those 
based on the hydrodynamic relationship among c, the 
spin stiffness, and the magnetic susceptibility BE may 
then be more convenient in practice. 


E. Width of the Classical Critical Region 

A key question raised by the experiments on TICUCI 3 
[5j concerns the width of the region close to Tjv where 
classical critical scaling applies. It was found that this 
width, W ~ 0.2Tjv, is essentially constant when normal¬ 
ized by Tn- We employ scaling arguments to show that 
the normalized width, W = W/Tn, should indeed be a 
constant with only a weak logarithmic correction in 3+1 
dimensions. 

For fixed temperature T, we consider the correlation 
length, defined in terms of the approach of the spin-spin 
correlation function to its asymptotic long-range value, 
?n 2 , when approaching the critical coupling ratio g{T) 
from the ordered side. This quantity has an initial diver¬ 
gence governed by the 4D QCP, 

£(T) = U(T) ~ \g c (T = 0) - g(T)]- v *, (34) 

with mean-field exponent iq = 1 / 2 , because the temporal 
thickness L T far exceeds £ and the system cannot sense 
its finite temporal extent, behaving as at T = 0. At the 
point where £ reaches L T , the behavior crosses over to a 
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3D scaling form, 

£{T) = UT) ~ [gJT) - g(T)}~^, (35) 

where ^3 « 0.70 is the 3D 0(3) exponent. Without 

logarithmic corrections, the temporal length is simply 
L t oc 1/T (more precisely, L T = L/c), but at the up¬ 
per critical dimension this relationship is modified by a 
logarithmic factor, 

L t ~ | ln(T)|®/T, (36) 

which is obtained by generalizing the classical result of 
Kenna [25). For the 4D 0(3) universality class, q = 1/4 
[ 2302 , and the correlation length itself also has a loga¬ 
rithmic correction, 

{g c -g)~ Ui \^{g c -g)V, (37) 

with v = 5/22. The quantum-classical crossover taking 
place when £ sa L T therefore corresponds to 

I ln(T)|®/T ~ \g c — g)]~ Vi \ hi(<7 c — g)]”, (38) 


which, by converting to a temperature-dependence and 
keeping only the leading logarithm, yields the crossover 
temperature 

T*(g)~(g c -g) 1/2 \Mg c -g)r°. (39) 

Using our result for T^(g) [Eq. ©] , the width of the 
classical critical region on the ordered side of the transi¬ 
tion is therefore 

W(g) = ~ l-a\Hg c -g)\^/\ (40) 

with a constant a, whose calculation requires further con¬ 
siderations, and a small exponent 


q-v 


T 


1/4-5/22 

3/11 


1 / 12 , 


(41) 


on the logarithm. This very weak dependence explains 
the near-constant behavior found for TICUCI3 ]§]. On 
general grounds we expect the width of the classical crit¬ 
ical regime on the other side of the transition to scale in 
the same way. 


VII. SUMMARY 

We have provided a direct and non-perturbative verifi¬ 
cation of the existence and nature of multiplicative loga¬ 
rithmic corrections to scaling at the quantum phase tran¬ 
sition for three-dimensional dimerized quantum Heisen¬ 
berg antiferromagnets. These systems correspond to the 
</> 4 field theory of an 0(3) quantum field in 3 + 1 dimen¬ 
sions, which is the upper critical dimension (D c = 4) for 
all models with O (TV) universality. With the exception 
of the Ising model (N = 1) |50| . no such demonstration 


exists to date, despite a significant body of analytical 
and numerical work on quantum criticality in dimerized 
quantum antiferromagnets. 

Our results are obtained from large-scale quantum 
Monte Carlo calculations based on state-of-the-art simu¬ 
lation techniques and detailed finite-size-scaling analysis. 
These enabled us to extract the precise logarithmic cor¬ 
rections to the leading critical properties at the quantum 
phase transition from a non-magnetic state of dominant 
dimer correlations to a Neel-ordered antiferromagnetic 
state. Specifically, we have obtained the multiplicative 
logarithmic corrections to the mean-field behavior of the 
order parameter, the zero-temperature staggered mag¬ 
netization (to s ), on the control parameter, the coupling 
ratio g. We have verified that these are governed by an 
exponent (3 = 3/11, a value we specify with numerical 
(statistical) precision under 3%, matching precisely the 
prediction of perturbative renormalization-group calcu¬ 
lations [211 1251 . 

No prediction was previously available for the analo¬ 
gous logarithmic correction to the Neel temperature, T/v- 
We have implemented a scaling Ansatz exploiting the 
known logarithmic corrections of other physical quanti¬ 
ties to obtain its form. Our prediction is that Tjv has 
exactly the same exponent in its logarithmic correction, 
r = 3/11, as the order parameter, and our numerical 
results for Tjsr(g) are in excellent agreement. We have 
thereby established an exact linearity between T/v and 
m s throughout the quantum critical regime. We have 
also demonstrated a different kind of logarithmic correc¬ 
tion, in the size-dependence of the staggered magnetic 
susceptibility at the four-dimensional quantum critical 
point, where we verify the predicted N- independent scal¬ 
ing form [35 ]. 

The numerical task of finding logarithmic corrections is 
not a straightforward one. We have established that the 
appropriate scaling regime is | g — g c \/gc 1$ 0.2. Within 
this region, obtaining reliable evidence for logarithmic 
corrections is critically dependent on having many high- 
precision data points at very small values of \g — g c \, 
which mandates accurate calculations at large system 
sizes. After establishing the location of the critical point 
to approximately one part in 10 5 , g c = 4.83704(6), we 
were able to obtain highly accurate extrapolations of the 
physical observables for coupling ratios as close to g c as 
4.834, i.e. with | g — g c \/g c — 0.0006. This required work¬ 
ing with linear system sizes as large as L = 48, mean¬ 
ing a system containing N = 2L 3 = 221184 interacting 
spins, and at temperatures as low as T = 1/(2 L) = 1/96. 
From this perspective, it becomes obvious why previ¬ 
ous studies m 121 ed, with only a handful of data 
points in the quantum critical regime (none closer than 
| g — g c \/g c = 0.02) were not able to find any meaningful 
evidence for logarithmic corrections. 

Our results are directly relevant to the pressure- 
induced quantum phase transition in TlCuCU [6HH]- De¬ 
tailed experiments on this material by elastic and inelas¬ 
tic neutron scattering have measured the staggered mag- 
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netization, the Neel temperature, the gap of the quan¬ 
tum disordered phase, and the magnetic excitation spec¬ 
trum on both sides of the transition. On the assump¬ 
tion that the leading dependence of the control param¬ 
eter (the ratio of antiferromagnetic superexchange pa¬ 
rameters) is linear in the applied pressure, both m s and 
Tpj show good mean-field exponents and a close linear 
relation over much of the accessible pressure range [5]. 
On the grounds that the available data follow mean-field 
scaling around the quantum critical point, it cannot be 
argued that they provide any evidence for logarithmic 
corrections, although the size of the experimental errors 
and the shortage of data very close to the QCP certainly 
mean they cannot be excluded. 

It has been argued very recently [53], based on a field- 
theoretic treatment, that the apparent suppression of m s 
and T/v visible in the experimental data for TICUCI3 
rather far from the QCP (at pressures 2-4 times the criti¬ 
cal pressure) arises due to logarithmic scaling of the cou¬ 
pling constant. In Ref. J5] it was assumed that these 
effects are in fact a consequence of departures from the 
quantum critical scaling regime, evident also in the vio¬ 
lation of linearity between m s and T/v beyond the point 
where the order parameter is 60% of the classical mo¬ 
ment. Although a direct comparison with our results is 
not possible without a microscopic treatment of the re¬ 
lationship between the applied pressure and the control 
parameter, a similar downturn is visible, and better de¬ 
scribed by including the multiplicative logarithmic cor¬ 
rections, beyond \g — g c \/gc ~ 0.05 in our Figs. [4] and [9] 
It is also tempting to relate the Tjv(m s ) curve of TICUCI3 
0to ourFig.^b), where the extended linear regime is 
followed by an upturn deep inside the Neel phase, which 
was interpreted m as the breakdown of the quantum- 
thermal decoupling (see below) due to a large density 
of thermally excited magnons when T/v is high. Finally, 
we have also provided a theoretical explanation for the 
shape of the classical critical scaling “fan” around T/v (p) 
observed in TICUCI3 [8], by showing that its width scales 
linearily with T/v, modified by a logarithmic correction 
with a very small exponent of 1/12, which would vary ex¬ 
tremely weakly over the experimental pressure window. 

Although many dimerized S = 1/2 systems with an¬ 
tiferromagnetic interactions are known, and many field- 
induced quantum phase transitions have been studied, 
few have yet been found to be close to quantum critical 
points at zero field under pressure. Our results shed light 
on the experimental challenges inherent in finding loga¬ 
rithmic corrections, but also provide evidence that their 
detection is actually possible. While important theoreti¬ 
cal questions remain to be addressed in lower dimensions, 
logarithmic corrections are of little relevance away from 
D c . Another challenge for both experiment and numeri¬ 
cal simulation would be to investigate the exponents and 
corrections for different N, meaning for systems of Ising 
and XY spins. A related experimental possibility would 
be to realize the N = 2 situation in a gas of ultra-cold 
bosons on an optical lattice. The unfrustrated dimer¬ 


ized antiferromagnet is a bipartite lattice and thus can 
be treated exactly as a system of hard-core dimer bosons, 
with the dimerized phase corresponding to the Mott in¬ 
sulator and the antiferromagnet to the superfluid (a state 
of long-range inter-site coherence); the symmetry broken 
is U(l), which is equivalent to XY. Although these ex¬ 
periments have not yet been realized in sufficiently large 
three-dimensional gases of cold bosons, the very fine pa¬ 
rameter control possible in cold-atom systems offers an¬ 
other candidate route for the experimental observation 
of logarithmic corrections to scaling. 

One of the points made in Ref. [8] was that, although 
quantum critical phenomena are universal, obeying scal¬ 
ing forms determined only by macroscopic properties of 
the system such as the dimensionality and the symmetry 
of the order parameter, their experimental observation 
depends crucially on non-universal prefactors. For quan¬ 
tum critical excitations, this is the ratio of the width of 
an excitation to its energy, and is a quantity determined 
entirely by microscopic details. For both static and dy¬ 
namic properties, the key figure of merit is the width of 
the quantum critical regime, and for this we have ob¬ 
tained a quantitative result not previously available by 
any other technique, \g — g c \/g c < 0.2. In as much as one 
may generalize from the dimensionality and geometry of 
the double cubic lattice, this 20% criterion dictates the 
necessary proximity to the quantum critical point for the 
observation of strict quantum critical scaling, including 
logarithmic corrections. 

Our demonstration of linearity between m s and T n 
in the (3 + 1)D Heisenberg antiferromagnet lies beyond 
any results previously predicted by analytical methods. 
What we have demonstrated explicitly for several quan¬ 
tities is the presence of expected logarithmic corrections, 
but their cancellation between m s and T/v was not antic¬ 
ipated. However, in parallel to our scaling argument for 
the logarithmic corrections to T/v, Scammell and Sushkov 
have recently arrived at the same conclusion from a dif¬ 
ferent starting point [53]. A key outstanding question is 
whether the linearity of T/v(m s ) is in fact a more funda¬ 
mental property of the system than arguments made at 
the semiclassical and mean-field levels suggest. Qualita¬ 
tively, the origin of linearity is thought m to lie in the 
effective decoupling of the classical and quantum fluctu¬ 
ations, which is applicable for all coupling ratios both 
outside and inside the QC regime. Its observation here 
implies the enduring independence of quantum and ther¬ 
mal fluctuations at the O(N) transition with D = D c 
for any N. Efforts to study the relationship between 
the T = 0 order parameter and the critical temperature 
in systems with different universality classes would shed 
light on this matter. 
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Appendix: Crossing-point scaling in the presence of 
logarithmic corrections 

According to the general hypothesis of finite-size scal¬ 
ing (FSS), verified by renomalization-group techniques 
(for a review see Ref. 52]), the dependence on system 
size of a physical quantity in the neighborhood of a crit¬ 
ical point can be described by the function 

Q(t, L) = + 0(L- U , r“% (A.l) 

where t is the distance to the critical point, i.e. t = 
\T — T c \ for a classical or t = \g — g c \ for a quantum phase 
transition, k is the critical exponent for the quantity in 
question in the thermodynamic limit, Q(t) ~ |f| _K , and 
the subleading exponent w originates from an irrelevant 
scaling field. In a fully rigorous treatment, L in f(£/L) 
should be replaced by £l(0 ), which is the finite-size cor¬ 
relation length at the critical point (t = 0 ) and thus the 
relevant length scale for FSS. In the absence of logarith¬ 
mic corrections, £l(0) ~ L. 


The leading term in Eq. (A.l) is the asymptotic FSS 


and the second term expresses the correction to scaling. 
For the Binder ratio, Q(g,L) = R 2 (g,L), which is a di¬ 
mensionless “invariant,” the asymptotic scaling has ex¬ 
ponent k = 0. Correction terms remain present, and at 
the crossing point t* for two system sizes L\ and L 2 one 
has 

R 2 (t*,L 1 ) = R 2 (t*,L 2 ). (A.2) 

Without logarithmic corrections, 

m/L) = h(tL V") (A.3) 

and hence 

R 2 {t,L) = a + btL 1/v + cL~ u + ... (A.4) 

The crossing point, t*, can be determined for (Li, L 2 ) as 
1 - s~ u 


t* 


s 1 /" - 1 


- —UJ—l/l' 
J 1 5 


(A.5) 


If the logarithmic correction to the correlation length 
is taken into consideration, 


£~f"|hitr 

and, according to Refs. ESI HZ], 

£ l ( 0 )~L 1 iPL 


(A. 6 ) 


(A.7) 


is n ow th e relevant FSS length scale. On substituting 
Eq. (A.71 into both the asympt otic F SS term f(£/L) and 
the subleading term L - ^, Eq. (A.5) becomes 


f* i* ~ s TT 7 Ary C _1/ » ( A - 8 ) 


where 


s = 


q 2 (0) = L 2 \n%L 2 ) 
^(0) Liln^LO’ 


(A.9) 


which also approaches a constant as Li,L 2 —> oo. 

There is no straigtforward inversion of Eq. ( |A .8 1 to 
obtain an exact expression for t*. However, for the lead¬ 
ing logarithmic correction it is sufficient to substitute 
Eq. (A.5) into the logarithmic part of (A. 8 ), which yields 


In t* ~ c + In ( - - ) — (w + 1/v) In L, (A.10) 


a quantity approximately proportional to In L when L is 
large (consider s = (L + 2)/L = 1 + 2/L). The le ading 
scaling behavior, obtained on replacing Inf* in Eq. (A. 8 ) 
by In L, is 


1 — UJ 

t* ~ - L -U-l/u ln -qu-q/v L ln O/u L 


S'l/v _ 1 

1 - s'- u 


L -ui-l/v ln c L 


s n / u - 1 

where the exponent c is given by 


„ v-q „ 
c =- quj. 


(A.ll) 


(A.12) 


Replacing t by q — q r and taking the large-L limit such 
that s' —^ 1, Eq. dATll yields 


9c 


(L) = g c + aL- (1/, ' +u>) In & L, 


(A.13) 


which is Eq. R 8 | in Sec. Ill A However, if there is no 


logarithmic correction to the subleading term L w , the 
second term in Eq. (A. 12) is absent, and 


where s = L 2 IL\ : and is either constant (L 2 = aL\ with 
a > 0 ) or approaches unity (L 2 = L\ + 2 ) as the system 
size goes to infinity. 



v 


(A.14) 
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